/* peq.reajs
 * 
 * Copyright (c) 2012, Michael Gruhn
 * All rights reserved.
 * 
 * Permission to use, copy, modify, and/or distribute this software for any
 * purpose with or without fee is hereby granted.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY, FITNESS AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
 * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
 * PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
 * THIS SOFTWARE.
 */

/* Parametric EQ after [1].
 */
desc: Parametric EQ

slider1:1000<20,20000,1>Frequency (Hz)
slider2:0<-12,12,1>Gain (dB)
slider3:1<0,4,0.1>Bandwidth (Oct)
slider6:1<0,1,{digital,analog}>Curve shape
slider7:0<0,1,{1/2 Gain, 3dB}>Width at

@slider
G0 = 1;
G = pow(10,slider2/20);

slider2 ? (
	GB = pow(10,slider2/2/20);
	slider7 ? (
		abs(slider2/2) > 3 ? (
			GB = pow(10,3*sign(slider2)/20);
		);
	);
	
	w0 = 2*$pi*slider1/srate;
	tmp = (slider1*2^slider3-slider1)/srate * 2;
	tmp = (exp(tmp)-exp(-tmp))/(exp(tmp)+exp(-tmp)) * 0.5 / 2;
	Dw = 2*$pi*tmp;
	
	F = abs(G^2 - GB^2);
	G00 = abs(G^2 - G0^2);
	F00 = abs(GB^2 - G0^2);
	num = G0^2 * (w0^2 - $pi^2)^2 + G^2 * F00 * $pi^2 * Dw^2 / F;
	den = (w0^2 - $pi^2)^2 + F00 * $pi^2 * Dw^2 / F;
	G1 = slider6 ? sqrt(num/den) : 1;
	G01 = abs(G^2 - G0*G1);
	
	G11 = abs(G^2 - G1^2);
	F01 = abs(GB^2 - G0*G1);
	F11 = abs(GB^2 - G1^2);
	W2 = sqrt(G11 / G00) * tan(w0/2)^2;
	DW_ = (1 + sqrt(F00 / F11) * W2) * tan(Dw/2);
	C = F11 * DW_^2 - 2 * W2 * (F01 - sqrt(F00 * F11));
	D = 2 * W2 * (G01 - sqrt(G00 * G11));
	A = sqrt((C + D) / F);
	B = sqrt((G^2 * C + GB^2 * D) / F);
	
	b0 = (G1 + G0*W2 + B) / (1 + W2 + A);
	b1 = -2*(G1 - G0*W2) / (1 + W2 + A);
	b2 = (G1 - B + G0*W2) / (1 + W2 + A);
	a0 = 1;
	a1 = -2*(1 - W2) / (1 + W2 + A);
	a2 = (1 + W2 - A) / (1 + W2 + A);
):(
	a0 = a1 = a2 = b1 = b2 = 0;
	b0 = 1;
);

redraw = 25;

@sample
x0 = spl0;
y0 = b0*x0 + b1*x1 + b2*x2 - a1*y1 - a2*y2;
y2 = y1; y1 = y0; x2 = x1; x1 = x0;
spl1 = spl0 = y0;

@gfx
(redraw+=1) > 25 ? (
redraw = 0;

gfx_r=gfx_g=gfx_b=0; gfx_a=1;
gfx_x=gfx_y=0;
gfx_rectto(gfx_w,gfx_h);

gfx_r=gfx_g=gfx_b=0.25; gfx_a=1;
gfx_y = 0; gfx_x = (log(20)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(30)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(40)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(50)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(60)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(70)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(80)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(90)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=gfx_g=gfx_b=0.5; gfx_a=1;
gfx_y = 0; gfx_x = (log(100)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=gfx_g=gfx_b=0.25; gfx_a=1;
gfx_y = 0; gfx_x = (log(200)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(300)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(400)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(500)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(600)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(700)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(800)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(900)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=gfx_g=gfx_b=0.5; gfx_a=1;
gfx_y = 0; gfx_x = (log(1000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_r=gfx_g=gfx_b=0.25; gfx_a=1;
gfx_y = 0; gfx_x = (log(2000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(3000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(4000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(5000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(6000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(7000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(8000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(9000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=gfx_g=gfx_b=0.5; gfx_a=1;
gfx_y = 0; gfx_x = (log(10000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=0.5; gfx_g=gfx_b=0; gfx_a=1;
gfx_y = 0; gfx_x = (log(srate/2)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=gfx_g=gfx_b=1; gfx_a=1;
gfx_x = 0; gfx_y = gfx_h - (0 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);

gfx_r=gfx_g=gfx_b=0.5; gfx_a=1;
gfx_x = 0; gfx_y = gfx_h - (6 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
gfx_x = 0; gfx_y = gfx_h - (-6 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);

gfx_r=gfx_g=gfx_b=0.25; gfx_a=1;
gfx_x = 0; gfx_y = gfx_h - (3 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
gfx_x = 0; gfx_y = gfx_h - (-3 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
gfx_x = 0; gfx_y = gfx_h - (9 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
gfx_x = 0; gfx_y = gfx_h - (-9 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);


g_x = 0;
gfx_x = 0;
gfx_y = gfx_h/2;
loop(gfx_w-1,
//	g_w = 2*$pi*(((g_x/gfx_w)^g_s*19980+20)/srate);
	tmp = 10^(g_x/gfx_w*3+log(20)/log(10));
	g_w = 2*$pi*tmp/srate;
	g_phi = 4*sin(g_w/2)^2;
	g_db = 10*log( (b0+b1+b2)^2 + ( b0*b2*g_phi - (b1*(b0+b2) + 4*b0*b2) )*g_phi )/log(10)
	- 10*log( (a0+a1+a2)^2 + ( a0*a2*g_phi - (a1*(a0+a2) + 4*a0*a2) )*g_phi )/log(10);
	
	g_y = gfx_h - (g_db + 12)/24 * gfx_h;
	
	gfx_r=0;gfx_g=1;gfx_b=0; gfx_a=1;
	
	gfx_lineto(g_x+1,g_y,0);
	
	gfx_x=g_x;
	gfx_y=g_y;
	g_x+=1;
);

);

/* 
 * References:
 * [1] Digital Parametric Equalizer Design With Prescribed Nyquist-Frequency Gain
 *     Sophocles J. Orfanidis
 */
